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FREQUENCY  TRACKING  WITH  HISTOGRAM-PMHT 


1.  INTRODUCTION 


Tracking  an  unknown,  nonstationary  signal  in  time  and  frequency  against  a  noisy,  non¬ 
stationary  background  on  an  intensity-modulated  sensor  display  is  a  difficult  problem.  Tradi¬ 
tional  techniques  involve  thresholding  the  sensor  data  and  treating  the  exceedences  as  point 
measurements  that  are  subsequently  fed  to  a  tracking  algorithm.  Choosing  this  threshold 
is  a  challenge  in  itself,  and  even  then  it  is  typically  subject  to  a  prescribed  probability  of 
detection  and  probability  of  false  alarm.  In  reference  1,  Streit  derives  a  new  tracking  algo¬ 
rithm  that  uses  all  of  the  senor  data,  and  thus  avoids  thresholding  entirely.  The  fundamental 
premise  of  this  tracking  algorithm  is  that  losses  due  to  thresholding  the  sensor  data  can  be 
eliminated  completely  if  all  of  the  sensor  data  are  used  by  the  tracking  algorithm. 

Section  2  describes  this  new  tracking  algorithm,  referred  to  herein  as  histogram  proba¬ 
bilistic  multi-hypothesis  tracking,  or  histogram-PMHT;  in  particular,  the  key  aspects  of  its 
derivation  are  discussed.  A  concise  statement  of  the  histogram-PMHT  algorithm  applied 
to  frequency  tracking  is  given  in  section  3.  Two  examples,  one  involving  a  simulated  linear 
chirp,  and  one  involving  a  bowhead  whale  call  recorded  at  sea,  are  presented  in  section  4.  A 
summary  is  given  in  section  5. 


2.  HISTOGRAM-PMHT 

2.1  SIGNAL  AND  MEASUREMENT  MODELS 

The  histogram-PMHT  tracking  algorithm  is  intrinsically  a  multi-signal  tracking  algo¬ 
rithm,  and  is  based  on  a  stochastic  model  of  the  signals  and  the  noise  background.  It 
assumes  that  the  signals  and  the  noise  background  are  described  by  a  discrete  mixture  of 
continuous  distributions  in  which  the  noise  background  and  each  signal  are  represented  by 
a  unique  component,  or  set  of  components,  in  the  mixture.  The  sensor  data  are  indirect 
manifestations  or  realizations  of  this  underlying  distribution.  Since  in  most  practical  appli¬ 
cations  the  sensor  display  has  a  fixed  resolution  with  a  finite  number  of  cells  (e.g.,  discrete 
Fourier  transform  (DFT)  bins  in  a  time-frequency  “waterfall”  display),  a  discrete  distribu¬ 
tion  is  required  to  model  the  sensor  data.  The  approach  taken  in  histogram-PMHT  is  first 
to  quantize  the  real-valued  sensor  data  into  a  “pseudo-histogram,”  and  then  to  use  a  multi¬ 
nomial  distribution  to  model  the  counts  in  the  histogram  cells.  The  cell-level  intensities  of 
the  sensor  data  are  directly  proportional  to  the  cell  counts  of  this  pseudo-histogram.  The 
goal  is  to  fit  the  underlying  mixture  distribution  to  this  pseudo-histogram  at  each  scan;  that 
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is,  to  estimate  the  parameters  of  the  mixture  distribution  that  maximize  the  likelihood  of 
the  pseudo-histogram  at  each  scan.  The  theoretical  framework  of  PMHT  (see  reference  2) 
is  used  to  “assign”  the  pseudo-histogram  samples  to  the  mixture  components,  and  to  link 
these  mixture  distributions  across  scans  with  a  dynamical  signal  model. 


2.2  EXPECTATION-MAXIMIZATION 

The  expectation-maximization  (EM)  method  formalized  by  Dempster,  Laird,  and  Rubin 
(see  reference  3)  is  a  powerful  method  for  estimating  the  parameters  of  mixture  distributions, 
and  is  the  method  used  to  solve  the  likelihood  equation  in  histogram-PMHT.  The  method  of 
EM  is  particularly  well  suited  to  so-called  “missing”  data  problems;  that  is,  problems  in  which 
the  parameters  of  interest  are  comparatively  straightforward  to  estimate  if  the  observed  data 
set  is  augmented  with  certain  unobserved  data.  The  basic  strategy  of  EM  is  to  include  the 
missing  data  as  random  variables  in  the  likelihood  function,  take  the  expectation  of  the  log- 
likelihood  with  respect  to  the  missing  data  conditioned  on  the  observed  data,  and  maximize 
the  resulting  expression,  termed  the  auxiliary  function.  The  EM  method  requires  both  the 
specification  of  the  incomplete  (observed)  data  likelihood  function,  and  the  complete  data 
(observed  data  plus  missing  data)  likelihood  function.  The  incomplete-  and  complete-data 
likelihood  functions  for  histogram-PMHT  are  described  in  the  next  two  sections. 


2.3  INCOMPLETE-DATA  LIKELIHOOD  FUNCTION 


Let  Zt  denote  the  sensor-data  measurement  vector  at  time  t, 

Zt  =  {zti,-  •  •  ,ztL},  t  =  l,...,T,  (2-1) 

where  z u  is  the  output  of  the  sensor  at  time  t  in  display  cell  £  (unaveraged,  short-time, 
magnitude-squared  Fourier  transform  data  versus  frequency  bins  in  this  application).  Let 
h2  >  0  be  a  specified  quantization  level,  and  let 


Nt  —  . . . ,  ti t  1, . . . ,  T, 

denote  the  quantized  measurement  vector  corresponding  to  Zt,  where 


is  the  greatest  integer  less  than  or  equal  to  Let 

L 

Nfz =y^,nu 

i=\ 


(2.2) 


(2.3) 


(2.4) 
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denote  the  total  cell  count  (sample  size)  at  time  t.  The  quantized  data  vector  Nt  is  used  as  an 
intermediate  variable  in  the  derivation  of  the  histogram-PMHT  algorithm;  at  an  appropriate 
point  in  the  derivation,  the  measurement  vector  Zt  is  recovered  in  the  limit  as  Ti2  — »  0. 


It  is  assumed  that  the  quantized  data  vector  Nt  has  a  multinomial  distribution  consisting 
of  Nt£  independent  draws  (with  replacement)  on  L  “categories”  with  probabilities 


where 


and 


SS 

II 

i— > 

N 

(2.5) 

Pi(Xt)  =  f f(r-,Xt)dr , 

(2.6) 

h 

L 

P(Xt)  =  Y^Pi(Xt). 

t= i 

(2.7) 

In  equation  (2.6),  f(r-,Xt )  is  the  underlying  mixture  density  of  the  signals  and  the  back¬ 
ground  noise,  and  Xt  denotes  the  parameter  vector  of  the  mixture  density  at  time  t  (minus 
the  mixing  proportions,  which  are  implicit  in  all  of  the  likelihood  expressions  in  the  sequel). 
The  assumption  of  a  multinomial  distribution  for  the  quantized  data  vector  Nt  implies  that 
the  counts  Nt  form  a  histogram  with  L  cells  and  a  sample  size  of  iVtS,  where  the  samples 
are  independent  and  identically  distributed  with  probability  density  function  f{r\Xt)  (see 
reference  4  for  more  on  the  multinomial  distribution). 


Let  N  =  {iV1, . . . ,  Nt}  denote  the  collection  of  quantized  measurement  vectors,  and  let 
X  =  denote  the  set  of  mixture  parameters  to  be  estimated.  Assuming  the 

scans  are  independent,  the  incomplete  data  likelihood  function  for  N  is  given  by 

T 

X)  =  JJ  Pi„c(N,-,  X, ),  (2.8) 

t=l 


where 


Pinc(Nt]  Xt) 


ATffi! 

nn\---ntL\ 


n 


P(X.) 


(2.9) 


In  reference  1,  a  Bayesian  model  for  the  mixture  parameters  X  is  adopted.  If  ps{X) 
denotes  the  a  priori  density  of  X ,  then  the  incomplete  data  likelihood  function  is  given  by 

Pinc(N,X )  =pE(X)pinc(N\X),  (2.10) 
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where  the  density  Pinc(N\X )  is  essentially  identical  to  the  density  (2.8),  and  differs  only 
in  its  statistical  interpretation  (Bayesian  conditioning  versus  parametric  dependence).  The 
derivation  of  the  prior  density  pz(X)  is  an  important  theoretical  development  in  histogram- 
PMHT.  In  short,  the  prior  needs  to  be  sufficiently  non-diffuse  so  that  the  synthetically 
generated  histogram  counts  N,  which  depend  on  the  arbitrary  quantization  level  fi2,  do  not 
overwhelm  the  prior  as  H 2  — >  0.  Under  the  usual  Markov  assumption  on  the  signal  states 
Xt,  Bayes  Theorem  gives 

T 

MX)  =S5,(X0)  (2-n) 

t= 1 

where  the  power  of  Ntz  is  necessary  to  account  for  the  artificial  abundance  of  quantized  data 
Nt  at  each  scan. 


2.4  COMPLETE-DATA  LIKELIHOOD  FUNCTION 

The  missing  data  in  histogram-PMHT  are  (1)  the  locations  of  the  samples  that  make  up 
the  pseudo-histogram,  and  (2)  their  mixture  component  assignments,  that  is,  the  components 
in  the  underlying  mixture  distribution  from  which  the  samples  are  drawn. 

Let  Cte  =  {Ctn,  •  •  • ,  C«nM}  denote  the  locations  of  the  samples  within  cell  i.  The  random 
variables  in  are  assumed  to  be  independent  and  identically  distributed  with  probability 
density  function  f(r\Xt)/Pe(Xt).  Furthermore,  let  Ct  =  (Cti,  •  •  • ,  Cti}  and  C  =  (Ci>  •  ■  • ,  Cr}- 

The  sample  density  is  assumed  to  be  the  mixture  density 

M 

1(t\X,)  =  Y,*tkGt(T\X.),  (2.12) 

fc= 0 


where  the  mixing  proportions  7 rt*  >  0  for  all  t  and  k ,  ^tk  =  1  for  all  t,  and  the  Gk{r\Xt) 

are  themselves  probability  density  functions.  For  component  irtoGo^Xt)-,  7rt0  represents  the 
fraction  of  the  total  power  due  to  the  background  noise,  and  Go(r\Xt)  models  the  cell-to-cell 
variation  of  the  background  noise.  Likewise,  for  components  7rtl  Gi(r\Xt), . . . ,  ntM  G& f(r\Xt), 
7 rtk  is  the  fraction  of  the  total  power  due  to  signal  k,  and  Gk(r\Xt)  models  the  variation  of 
signal  k  from  cell  to  cell.  The  mixture  model  (2.12)  assumes  that  the  signal  power  levels 
may  be  spread  across  more  than  one  cell  of  the  sensor  display. 

Let  Ku  =  {kui,  •  ■  ■ ,  kttnu}  denote  the  components  of  the  mixture  that  generated  the 
missing  variables  Cte  =  {C«i>  •  •  •  It  is  assumed  that  the  random  variables  in  Ku 

are  independent  and  identically  distributed.  Furthermore,  let  Kt  =  {Kt\, . . . ,  Kti}  and 
K  =  {Ku...,Kt}. 
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Including  the  missing  sample  locations  and  their  mixture  component  assignments,  the 
complete  data  likelihood  function  is  given  by 


T  L  nti 

Pcom(N,<;,K,X)=pE(X)  nnn  fkttr(Ct£r\Xt ),  (2-13) 

t=  1  1=1  r=l 

where 


fk(T\Xt)=irtkGk(r\Xty,  (2.14) 

that  is,  the  complete  data  likelihood  function  is  the  product  of  all  the  histogram  sample  den¬ 
sities  across  all  L  cells  and  all  T  scans,  scaled  by  the  prior  for  the  signal  mixture  parameters. 


2.5  E-STEP 


The  auxiliary  function  of  the  EM  method,  denoted  here  by  \&,  is  defined  as  the  conditional 
expectation  with  respect  to  the  missing  data  of  the  logarithm  of  the  complete  data  likelihood 
function  given  the  observed  data  and  the  current  values  of  the  parameters  X,  denoted  X': 

*(X,X')=EiK[\ogPcom(N,C,K,X)\N,X'],  (2.15) 

where  E^K  denotes  expectation  with  respect  to  the  missing  data.  The  mechanics  of  the 
E-step  for  histogram-PMHT  are  tedious  but  straightforward,  and  are  well  documented  in 
reference  1.  The  final  result  in  terms  of  the  pseudo-histogram  counts  N  is 


*(*.*')  =  10gPH.(X„)  + 


T  L  M 


******  A 

EEE  pTFT  /  'oeAM**)  (2-16) 


It  is  easily  shown  that  by  taking  the  limit 


¥  =  lim  h2  , 
h2-+  o 


(2.17) 


and  using  definition  (2.3)  for  the  quantization,  the  EM  auxiliary  function  ^  can  be  replaced 
by  a  function  of  the  unquantized  sensor  data  Z  =  {Zi, . . . ,  Zt }: 


nx,x')=j2wxk 

t=i 

T  L  M  , 

+  EEE-57VT  MtW)  log  Mt\X,)  dr,  (2.18) 

i=i  /=i  k=o 

where  ||Zt||  =  J2t=i  zti  is  the  -^i-norm  of  Zt. 
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2.6  M-STEP 


The  objective  of  the  M-step  is  to  maximize  the  auxiliary  function  ^  with  respect  to 
the  unknown  signal  parameters  X.  To  proceed,  application-specific  terms  in  the  auxiliary 
function  (2.18)  must  be  defined. 

Let 


Xt  —  3>tii  •  •  •  ■>  >  t  —  I?  -  *  *  >  T,  (2.19) 

where  Xto  is  the  background  noise  parameter  at  time  t,  and  Xtk  is  the  parameter  of  signal  k 
at  time  t.  Assuming  the  signals  are  independent  at  all  times, 

M 

k=0 

It  is  readily  shown  that  can  be  separated  into  two  terms,  one  involving  the  unknown 
mixing  proportions  7r  =  {w**},  and  one  involving  the  unknown  signal  parameters  X, 

T  M 

(2.21) 

t=l  k= 0 


where 


M 


pTF)  / G^r\x'tk)  dr 
k=0  U=1  Pt  ™ 


n'tk  logTTtifc, 


(2.22) 


and 

® kx 


+  £  £  Ife  f  ft(rfe)  logG(r|*»)  dr.  (2.23) 

The  updated  mixing  proportions  7rt  are  easily  obtained  at  each  time  t  by  maximizing  with 
respect  to  7rt  the  Lagrangian  equation  involving  \FtJr  and  the  constraint  7rt0-l-7rtiH f-7r tM  =  1- 

For  this  application,  linear  Gauss-Markov  processes  are  assumed  for  the  signals,  so  that 
for  k  =  1, . . . ,  M  and  t  =  1, . . . ,  T,  the  signal  process  models  are  given  by 


l,k)  —  ■M’iztk'i  F t— i,fc  Xt— i,fc>  Qt—  l,fc)» 


(2.24) 
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where  M{r\  C)  denotes  the  multivariate  normal  probability  density  function  with  mean 
H  and  covariance  matrix  C,  the  Ft-ilk  are  known  state  transition  matrices,  and  the  Qtk  are 
known  process  covariance  matrices. 

Additionally,  it  is  assumed  that  the  signal  components  in  the  mixture  distribution  are 
also  Gaussian,  and  that  the  means  of  these  Gaussians  are  linearly  related  to  the  states  of 
the  signals  k  =  1, . . . ,  M  at  times  t  =  1, . . .  ,T,  so  that 

Gk(r\xtk)  =  Af(r;  Htk  Xtk,  Rtk ),  (2.25) 

where  the  Htk  are  known  measurement  matrices,  and  the  Rtk  are  known  measurement  co- 
variance  matrices. 

Finally,  the  background  noise  distribution  in  this  application  is  assumed  to  be  uniform 
and  known  for  all  t,  so  that  the  Go(r|a:to)  terms  are  all  constants. 

With  these  assumptions,  it  is  shown  in  reference  1  that  for  X(k )  =  {xok,Xik, . .  .,xrk}, 
the  value  of  X  ( k )  that  maximizes  the  auxiliary  function  $>kx  for  each  signal  k  is  the  solution 
to  a  symmetric,  block-tridiagonal,  linear  system  of  equations,  and  that  this  system  is  most 
efficiently  solved  by  a  recursive  Kalman  smoothing  filter.  The  details  of  this  result  are 
omitted  here,  but  the  filter  steps  are  listed  explicitly  in  the  next  section. 


3.  FREQUENCY  TRACKER  ALGORITHM  STATEMENT 


In  the  case  of  frequency  tracking,  the  signal  parameters  of  interest  are  typically  instan¬ 
taneous  frequency  7*  and  instantaneous  frequency-rate  7t  at  time  t,  so  that  for  signal  k , 


Xtk  — 


(3.2) 


For  this  two-state  linear  Markov  model,  the  state  transition  matrices  Ft-i<k  and  the  process 
covariance  matrices  Qt-i,k  have  simple  forms: 


Ft-i,k  = 


1  A*-! 
0  1 


where  A*  is  the  elapsed  time  between  time  i  and  time  t  —  1,  and 

L  2^t-l  At_i 

where  the  qt-i,k  are  scale  factors  (see  reference  5). 


Qt-i,k  — 


(3.3) 


(3.4) 
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For  a  two-dimensional  time-frequency  display,  the  measurement  matrices  Htk  and  the 
measurement  covariance  matrices  Rtk  also  have  simple  forms: 

Htk  =  [1  0]  ,  (3.5) 

and 

Rtk  =  Ptk,  (3-6) 

where  Ptk  is  the  variance  of  signal  k  at  time  t. 

Let  {7^ }  be  the  set  of  estimated  mixing  proportions  and  {xj^}  the  set  of  estimated 
signal  states  at  the  i-th  EM  iteration.  At  the  beginning  of  the  algorithm  (the  0-th  iteration), 

the  mixing  proportions  {7^}  are  initialized  so  that  >  0  and  7 +  7rj^  H - h  7rj^  =  1 

for  t  =  1, . . . ,  T.  The  signal  state  sequences  {x^ ,  x^ , . . . ,  xj^}  are  initialized  with  nominal 
values  for  t  =  1,...,T.  For  iterations  i  =  0,1,2,...,  the  following  nine  quantities  are 
computed: 


1.  Component  bin  probabilities  for  t  =  1, . . . ,  T,  i  =  1, . . . ,  L,  and  k  =  0, 1, . 

0, 


,M: 


p< i+l)  _  Jl/£>  k 

*  =  1, 


(3.7) 


2.  Total  bin  probabilities  for  t  =  1, . . . ,  T  and  l  =  1, . . . ,  L: 

M 

pO+1)  _  V'  -.<0  P(i+1) 

ru  —  2-*/ 71  tk  ' 

k=0 


3.  Total  scan  probabilities  for  t  =  1, . . . ,  T: 

Pt(<+1)  =  pt{i+1) 

£=l 


(3.8) 


(3.9) 


4.  Bin  centroids  for  t  =  1, . . . ,  T,  £  =  1, . . . ,  L,  and  k  —  1, . . . ,  M: 


(3.10) 


5.  Synthetic  measurements  for  t  =  1, . . . ,  T  and  k  =  1, . . . ,  M: 


=(i+i)  _ 

ztk  ~ 


Etl 

[ztt 

(7&+1)/-P£+,)) 

z 

(i+l) 

tkl 

Zl. 

=1 

(3.11) 
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6.  Synthetic  measurement  covariance  matrices  for  t  =  1, . . . ,  T  and  k  =  1, . . . ,  M: 


6(*+i)  _ 

Ktk  — 


Rtk 


7T. 


tk 


V' L  (i+l)  (  p(*+l)  /p(»+l)\ 
Z-,£=lZt£  \rtk(  !  rtl  J 


(3.12) 


7.  Synthetic  process  covariance  matrices  for  t  =  0, 1, . . . ,  T  —  1  and  k  =  1, . . . ,  M : 

p(*+i) 


(3.13) 


8.  Estimated  mixing  proportions  for  t  =  1, . . . ,  T  and  k  =  0, 1, . . . ,  M: 

V r(*+1)  (  p(*+1)  /p(*+1)'\ 

l^i=i  zti  yrtu  !rti  J 


7r(*) 

_(*+!)  _  tk 
‘'tk  ~ 


\~\M  («)  /  p(i+l)  /p(i+l)\ 

A^k'=l  7rtfc/  2-*l=\  \rtk‘l  /  rtl  J 


(3.14) 


9.  Estimated  signal  states  for  t  =  0, 1, . . . ,  T  and  k  =  1, . . . ,  M,  using  (for  computational 
efficiency)  a  recursive  Kalman  smoothing  filter,  which  comprises  a  forward  filter  ini¬ 
tialized  at  time  t  =  0  with  y^l\k)  =  x^J  and  a  large  (diffuse)  state  covariance  matrix 

P^\k),  and  given,  for  t  =  0, 1, . . .  ,T  —  1,  by  the  recursions 


$(*)  =  FttP<;+l)(*)Fi  +  Qtt, 


(3.15) 

(3.16) 

(3.17) 

(3.18) 


and  a  backward  filter  initialized  at  time  t  =  T  with  Xj (&)  and  given,  for 
t  =  T  —  1,  T  —  2, . . . ,  1,  by  the  recursion 


and,  for  t  =  0,  by 


x0  k 


T(*+l) 

xlk  i 


(3.19) 


(3.20) 


where  the  asterisk  denotes  matrix  transposition. 


The  most  common  convergence  tests  for  termination  of  the  algorithm  are  based  on  the 
rate  of  increase  of  the  incomplete  data  likelihood  function.  Other  tests  are  based  on  the  rate 
of  change  of  the  estimated  parameters.  In  practice,  usually  some  combination  of  these  two 
tests  is  used. 
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4.  EXAMPLES 


4.1  LINEAR  CHIRP 


In  this  section,  a  two-state  histogram-PMHT  frequency  tracker  is  used  to  track  the  in¬ 
stantaneous  frequency  and  frequency-rate  of  a  low-frequency  linear  chirp. 

Consider  the  continuous-time,  complex-sinusoidal  signal  x(t )  with  constant  amplitude  A , 
constant  phase  4 i>,  and  time- varying  instantaneous  frequency  7 (t), 

x(t )  =  A  exp[z  27T7(t)  t  +  i<f>],  (4.2) 

where 


7(t)  =  7o  +  7o  t,  (4.3) 

70  is  the  nominal  frequency  in  Hz,  and  70  is  the  chirp-rate  in  Hz/s.  Converting  to  radians, 
equations  (4.2)  and  (4.3)  are  rewritten  as 

x(t)  =  A  exp[i  f2(t)  t  +  i  4>\,  (4.4) 

and 

fi(t)  =  Ho  +  flo  t ,  (4.5) 

where  fl(t)  is  the  instantaneous  frequency  in  rad/s,  fio  is  the  nominal  frequency  in  rad/s, 
and  f2o  is  the  chirp-rate  in  rad/s2. 


Without  loss  of  generality,  it  is  assumed  that  the  signal  is  observed  starting  at  time  t  =  0, 
and  that  estimates  of  the  instantaneous  frequency  and  frequency-rate  of  x(t)  are  required 
every  T  seconds.  To  avoid  aliasing,  the  signal  is  sampled  well  above  the  Nyquist  rate  for  the 
observation  period  of  interest  ST,  where  S  is  the  number  of  scans;  that  is 

n,  =  -^»2  n(ST),  (4.6) 

Is 


where  Cls  is  the  sampling  rate  in  rad/s,  and  Ts  is  the  sampling  period  in  seconds.  Using  n  to 
denote  the  sampling  index,  the  discrete-time  versions  of  equations  (4.4)  and  (4.5)  are  given 
by 


and 


x(nTs) 

x[n] 


A  exp[i  Q(riTs)  nTa  +  i  4>] 
A exp(i u[n\  n  +  i<f>), 


(4.7) 


Ts£l{nTa)  =  TsSlo  +  T,&onTs 
<j[n]  =  uo+uon, 
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(4.8) 


for  n  =  0, . . .  ,N  —  1,  where  u[n]  is  the  instantaneous  frequency  in  rad/sample,  u>o  is  the 
nominal  frequency  in  rad/sample,  and  uo  is  the  chirp-rate  in  rad/sample2  (see  reference  6 
for  details  on  time-sampling). 

For  simplicity,  it  is  assumed  that  the  signal  x(t )  is  corrupted  by  additive,  complex,  zero- 
mean  white  Gaussian  noise  with  variance  a2,  denoted  v(t)  ~  CA/"(0,  cr2);  that  is,  the  observed 
time-series  is  given  by 


y[n]  =  x[n)  +  u[n],  (4.9) 

where 

P{u[n]}  =  0,  (4.10) 

E{v[n]v[n  +  k]}  =  o2  5{k\,  (4.11) 

for  all  n  and  k.  An  observation  interval  of  length  T  seconds  and  a  sampling  period  of  length 
Ts  seconds  yields  a  time-series  of  length  T/Ts  +  1  samples  every  T  seconds,  for  a  total  of 
N  =  S(T/Ts  +  1)  samples. 

For  each  observation  interval,  the  signal-to-noise  ratio  (SNR)  77  is  defined  as  the  ratio  of 
the  average  signal  power  P  to  the  total  noise  power  a2, 

03  £  (4.12) 

For  the  signal  of  interest,  the  average  signal  power  P  =  A2,  so  77  =  A2 /a2. 

In  these  simulation  examples,  a  low-frequency  linear  chirp  with  amplitude  A  =  1,  phase 
<p  =  0,  nominal  frequency  70  =  20  Hz,  and  chirp-rate  70  =  0.125  Hz/s  is  sampled  every 
Ts  =  0.0125  second,  and  its  short-time  Fourier  transform  (STFT)  is  taken  every  T  =  1 
second  for  S  =  120  consecutive,  nonoverlapping  (unaveraged)  blocks  of  data  (scans)  to 
yield  a  spectrogram  of  length  ST  =  2  minutes.  This  sampling  period  corresponds  to  a 
sampling  rate  of  js  =  l/Ts  =  80  Hz,  which  is  well  above  the  Nyquist  rate  of  2  j(ST)  = 
2  [70  +  70 (ST)]  =  70  Hz  for  the  2-minute  observation  period.  Each  1-second  observation 
interval  yields  T / Ts  +  1  =  81  samples,  for  a  total  of  N  =  9720  samples  over  the  whole  data 
set. 


The  spectrogram  of  the  noise-free  (cr  =  0)  signal  is  shown  in  figure  1.  The  spectrograms 
for  signals  with  a  =  A  (0  dB),  3 A  (-9.5424  dB),  5 A  (-13.9794  dB),  and  6 A  (-15.5630  dB)  are 
shown  in  figures  3,  5,  7,  and  9,  respectively. 

The  STFTs  in  these  examples  are  computed  using  a  128-point  fast  Fourier  transform 
(FFT)  and  an  81-point  Hanning  window.  The  window  sequence  is  normalized  such  that  it 
sums  to  one;  this  scaling  allows  the  magnitude-squared  FFT  output  to  be  directly  interpreted 
as  power.  Since  for  a  P-point  DFT  the  spacing  between  DFT  frequencies  is  27r/P,  and  the 
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relationship  between  the  discrete-time  frequency  u  and  the  continuous-time  frequency  0  is 
u  =  flTs,  where  Ts  is  the  sampling  period,  the  frequency  resolution  is  given  by  2n/PTs  rad/s, 
so  the  DFT  bin-width  for  this  example  is  1/(128  x  0.0125)  =  0.625  Hz.  Only  the  positive 
half  of  the  spectrum  is  shown  in  figure  1,  and  only  the  first  L  =  64  bins  of  sensor  data  are 
considered  in  the  sequel. 

The  effect  of  the  window  is  to  smear  the  signal  in  time  and  frequency,  and  hence  add  to 
the  apparent  bandwidth  of  the  signal  induced  by  the  spectral  sampling.  While  windowing 
decreases  the  resolution  of  the  spectrogram,  it  would  seem  to  improve  the  fit  of  the  mixture 
model  to  the  data  in  this  case — the  broadening  effect  of  the  window  adds  Variance”  to  the 
signal.  Although  this  smearing  of  the  signal  by  the  window  is  a  purely  deterministic  effect, 
it  is  treated  as  a  stochastic  effect  by  the  histogram-PMHT  tracker;  that  is,  the  width  of  the 
signal  is  modeled  by  the  variance  of  the  signal  component  of  the  mixture. 

The  histogram-PMHT  frequency  tracker  outlined  in  section  3  was  implemented  as  a 
single-scan  (T  =  1),  single-signal  (M  =  1)  sequential  filter  and  was  applied  to  the  data  in 
figures  1,  3,  5,  7,  and  9,  on  a  per-scan  basis;  that  is,  the  scans  were  processed  sequentially, 
where  the  estimate  for  the  previous  scan  was  used  to  initialize  the  estimate  for  the  current 
scan  (see  reference  1).  The  values  used  for  qt-i  =  q  and  pt  =  p  in  equations  (3.4)  and  (3.6) 
are  listed  in  the  figure  captions. 

The  estimated  instantaneous  frequencies  are  shown  connected  by  a  solid  line  in  each 
figure.  The  tracks  become  increasingly  jagged  as  the  SNR  drops  and  the  signal  power  varies 
more  from  cell  to  cell.  The  jaggedness  of  the  tracks  is  also  a  function  of  the  process  noise:  the 
stochastic  signal  model  allows  for  further  excursions  from  the  deterministic,  constant-rate 
model  as  the  process  noise  scale  factor  q  is  increased,  resulting  in  larger  random  accelerations. 
There  is  an  intimate  relationship  between  the  values  of  the  filter  parameters  p  and  q ,  the  SNR, 
the  track  initialization,  and  the  tracker  behavior.  High  SNR  signals  will  tend  to  draw  the 
signal  component  of  the  mixture  closer  if  the  signal  is  within  the  effective  “gate”  determined 
by  p.  The  amount  by  which  the  signal  component  will  move  to  data  outside  the  range  of 
the  deterministic  process  model  depends  on  the  SNR,  the  value  of  q ,  and  the  track’s  local 
estimate  of  its  own  quality  (i.e.,  the  size  of  the  state  covariance  Pt\t  on  xt).  In  low  SNR,  the 
“inertia”  of  the  tracker  will  allow  it  to  coast  according  to  the  deterministic  process  model 
until  there  is  good  reason  (i.e.,  supporting  data)  for  it  to  change  course. 

The  track  in  figure  9  was  terminated  prematurely  because  it  had  lost  track.  Because  the 
nearest  peak  in  the  second  scan  of  data  is  far  to  the  left  of  where  the  track  was  initialized 
(at  21.875  Hz,  the  right  edge  of  bin  35),  a  large  innovation  was  generated  in  the  filter,  which 
led  to  a  large  negative  frequency-rate  estimate  from  which  the  tracker  never  recovered.  In 
the  absence  of  consistently  strong  data,  the  tracker  coasted  according  to  its  current  rate 
estimate,  with  minor  course  corrections  in  the  vicinity  of  strong  local  peaks.  This  example 
is  a  good  example  of  where  a  multi-scan  batch  filter  would  be  beneficial.  A  batch  of,  say, 
10  scans  in  this  example  would  probably  provide  enough  smoothing  for  an  acceptable  initial 
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track  estimate.  The  batch  estimate  could  then  be  used  to  either  initialize  the  next  batch, 
slid  forward  by  one  scan,  or  to  initialize  a  single-scan  sequential  tracker  like  the  one  used  in 
these  examples. 

Plots  of  the  instantaneous  frequency  and  frequency- rate  estimates  for  these  examples, 
along  with  the  true  frequencies  and  frequency-rate,  are  shown  in  figures  2,  4,  6,  and  8. 
The  apparent  small  (less  than  one  bin-width)  bias  in  the  frequency  estimates  is  thought 
to  be  due  more  to  a  combination  of  the  resolution  of  the  STFT,  and  the  known  bias  of 
the  windowed  periodogram  (see  reference  6),  than  to  the  tracker.  The  histogram-PMHT 
tracker  is  a  maximum  likelihood  estimator,  and  should  therefore  be  asymptotically  unbiased, 
efficient,  and  normal.  A  detailed  statistical  analysis  of  histogram-PMHT  performance  has 
not  yet  been  done,  and  is  left  for  future  work. 
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Scan 


Q  =  1  e-05  [  1/3  (I)3, 1/2  (I)2  ;  1/2  (1  )2,  4  ] ;  r  =  0.625  Hz  ;  PQ  =  [  (40  Hz)2,  0  ;  0,  (20  Hz/s)2  ] 


Figure  1.  Magnitude- Squared  STFT  of  Noise-Free  Linear  Chirp  (rj  =  oo  dB, 
p=  1  Bin  (0.625  Hz),  q  =  lx  10-5^  and  Estimated  Track 
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Frequency  (Hz)  Frequency  Rate  (Hz/s) 

Figure  2.  Estimated  and  True  Instantaneous  Frequencies  and  Frequency- Rates 
for  a  Noise-Free  Linear  Chirp  (r\  =  oc  dB,  p  =  1  Bin  (0.625  Biz),  q  =  lx  10-5^ 
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Q  =  le-05  [  1/3  (I)3, 1/2  (I)2  ;  1/2  (I)2, 4  ] ;  r  =  0.625  Hz ;  P0  =  [  (40  Hz)2,  0  ;  0,  (20  Hz/s)2  ] 


Figure  3.  Magnitude- Squared  STFT  of  Linear  Chirp  in  White  Gaussian  Noise 
(rj  =  0  dB,  p  =  1  Bin  (0.625  Hz),  q  =  lx  10~5j  and  Estimated  Track 
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RMS  =  0.41947Hz  (Scans  20  —  120) 


RMS  =  0.001 81 83Hz/s  (Scans  20  —  120) 


Figure  4.  Estimated  and  True  Instantaneous  Frequencies  and  Frequency- Rates 
for  a  Linear  Chirp  in  White  Gaussian  Noise  (rj  =  0  dB,  p  =  1  Bin  (0.625  Hz), 
q  =  lx  10"5; 
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Q  =  1e-05  [  1/3  (I}3, 1/2  (I)2 ;  t/2  (if,  4 ] ;  r  =  0.625  Hz  ;  PQ  =  [  (40  Hz)2, 0 ;  0,  (20  Hz/s)2 ] 


Figure  5.  Magnitude- Squared  STFT  of  Linear  Chirp  in  White  Gaussian  Noise 
(v  =  —9.5424  dB,  p  =  1  Bin  (0.625  Hz),  q  =  lx  10-5^  and  Estimated  Track 
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120 


RMS  =  0.003854Hz/s  (Scans  20  —  120) 


RMS  =  0.391 66Hz  (Scans  20  — 120) 


Figure  6.  Estimated  and  True  Instantaneous  Frequencies  and  Frequency-Rates 
for  a  Linear  Chirp  in  White  Gaussian  Noise  (rj  =  -9.5424  dB,  p=  1  Bin  (0.625 
Hz),  q=  lx  10"5; 
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Q  =  0.1  1 1/3  (I)3, 1/2  {1  )2 ;  1/2  (I)2,  4  ] ;  r  =  0.625  Hz  ;  PQ  =  [  (40  Hz)2,  0  ;  0,  (20  Hz/s)2 ] 


Bin 

Figure  7.  Magnitude- Squared  STFT  of  Linear  Chirp  in  White  Gaussian  Noise 
(rj  =  —13.9794  dB,  p  =  1  Bin  (0.625  Hz),  q  =  1  x  10-3^  and  Estimated  Track 
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RMS  =  0.55533Hz  {Scans  20  —  120) 


RMS  =  0.01 9242Hz/s  (Scans  20  —  1 20) 


Figure  8.  Estimated  and  True  Instantaneous  Frequencies  and  Frequency- Rates 
for  a  Linear  Chirp  in  White  Gaussian  Noise  (r)  =  —13.9794  dB,  p  =  1  Bin  (0.625 
Hz),  a  =  1  x  II 


Q  =  0.1  [  1/3  (I)3, 1/2  (I)2 ;  1/2  (if,  4  ] ;  r  =  0.625  Hz ;  P0  =  [  (40  Hz)2,  0  ;  0,  (20  Hz/sf  ] 


Figure  9.  Magnitude-Squared.  STFT  of  Linear  Chirp  in  White  Gaussian  Noise 
(t]  =  —15.5630  dB,  p  =  1  Bin  (0.625  Hz),  q  =  lx  10-3^  and  Estimated  Track 
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4.2  BOWHEAD  WHALE  CALL 


The  spectrogram  of  a  1.2-second  bowhead  whale  call  recorded  at  sea  is  shown  in  figure 
10.  This  signal  was  sampled  at  2500  Hz  and  processed  with  a  421-point  Chebyshev  window 
and  a  1024-point  FFT  with  a  bin  resolution  of  2.4414  Hz.  Values  of  p  =  3  bins  (7.3242  Hz) 
and  q  =  lx  10“3  were  used  to  track  this  signal.  The  estimated  instantaneous  frequencies 
are  shown  connected  by  a  solid  line  in  figure  10. 

The  intent  of  this  example  is  to  show  the  ability  of  the  tracker  to  track  very  complex,  non¬ 
stationary  signals,  and  the  ability  of  the  stochastic  process  model  to  accommodate  dynamics 
of  higher  order  than  the  deterministic  process  model.  The  benefit  of  including  process  noise 
is  the  ability  to  track  high  dynamics  with  a  reduced  parameter  set.  However,  for  too  low  a 
model  order,  the  burden  of  tracking  the  signal  dynamics  falls  squarely  on  the  shoulders  of 
the  process  noise,  and  the  result  for  high  process  noise  may  be  very  jagged  tracks  or,  for  low 
processes  noise,  tracks  whose  dynamics  lag  those  of  the  actual  track.  The  appropriate  model 
order  and  process  noise  level  are  of  course  application  dependent. 


5.  CONCLUDING  REMARKS 


The  histogram-PMHT  tracking  methodology  avoids  the  thorny  issues  of  thresholding  the 
sensor  data  to  provide  measurements  for  a  point  tracker  by  (1)  using  all  of  the  sensor  data, 
and  (2)  modeling  the  signals  as  potentially  distributed  over  several  cells.  The  histogram- 
PMHT  signal  model  is  a  stochastic  model  in  that  the  signal  centers  are  modeled  as  the 
component  means  of  a  discrete  mixture  distribution;  the  signal  “bandwidths”  are  modeled 
by  the  variances  of  the  mixture  components. 

The  interpretation  of  the  sensor  output  as  a  pseudo-histogram  plays  an  important  role  in 
the  derivation  of  histogram-PMHT.  The  real-valued  cell  outputs  are  quantized  and  treated 
as  the  cell  counts  of  a  pseudo-histogram  whose  distribution  is  multinomial;  the  underlying 
density  of  this  multinomial  distribution  is  the  mixture  density  of  the  signals  and  the  back¬ 
ground  noise.  It  is  a  remarkable  fact  that  the  sensor  data  are  recovered  in  the  algorithm  as 
the  quantization  level  h 2  is  taken  to  zero. 

Several  extensions  to  histogram-PMHT  exist,  and  some  of  these  are  developed  in  refer¬ 
ence  1.  For  instance,  the  measurement  covariance  matrices  Rtk  and  the  process  covariance 
matrices  Qtk  can  be  treated  as  unknowns  and  estimated  from  the  data.  Also,  some  of  the 
unknowns  can  be  linked  together  parametrically  or  held  fixed,  if  either  option  makes  sense 
in  the  application.  For  example,  the  process  covariance  matrices  Qtk  and  the  mixing  propor¬ 
tions  7 Ttk  can  be  held  constant  for  each  signal  over  all  scans  such  that  Qtk  =  Qk  and  tt**  =  i r* 
for  all  t  and  k . 
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Q  =  0.001  [  1/3  (I)3, 1/2  (I)2 ;  1/2  (if,  4  ) ;  r=  7.3242  Hz ;  P0  =  [  (40  Hz)2,  0 ;  0,  (20  H z&f  ] 


Figure  10.  Spectrogram  of  Bowhead,  Whale  Call 
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Another  important  theoretical  development  in  the  derivation  of  histogram-PMHT  not 
discussed  in  this  paper  is  the  use  of  a  negative  multinomial  distribution  (see  reference  4)  to 
model  sensor  cells  in  which  no  data  are  collected.  The  negative  multinomial  functions  as  an 
interpolator  in  this  case,  serving  to  restore  the  missing  cell  counts  in  the  pseudo-histogram. 
In  this  capacity,  the  negative  multinomial  may  be  useful  to  reduce  “edge  effects”  that  may 
bias  estimates  when  the  tails  of  signal  components  extend  beyond  the  ends  of  the  sensor 
display.  The  reader  is  referred  to  reference  1  for  further  details  on  the  use  of  the  negative 
multinomial  in  histogram-PMHT. 
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